#include <stdio.h>
#include <math.h>
#include <gmp.h>
#include <sys/time.h>
#include <time.h>

struct timeval inicio, final;

int main() {
	gettimeofday(&inicio, NULL);
	
	int i = 0, time;
	// Fator 2^(2*k+3) a ser multiplicado do calculo de an, k=0 resulta em 8 na primeira iteracao
	unsigned long int twoexp = 8;
	FILE *file;
	// O arquivo se chama "Borwein" e sera criado na mesma pasta do .c
    file = fopen("Borwein", "w+");
		
	// Precisao de 10m casas decimais, 10m * log2(10)
	mpf_set_default_prec(33219280);
	
	// Declara e inicializa variaveis
	mpf_t a, y, an, yn, temp;
	mpf_init(a);
	mpf_init(y);
	mpf_init(an);
	mpf_init(yn);
	mpf_init(temp);
	
	// Inicializa a0 e y0
	mpf_sqrt_ui(y, 2);
	mpf_mul_ui(a, y, 4);
	mpf_sub_ui(y, y, 1);
	mpf_ui_sub(a, 6, a);
	
	while(1) {
		printf("Iter n. %d\n", ++i);

		// yn = (1-(1-y^4)^(1/4))/(1+(1-y^4)^(1/4))
		mpf_pow_ui(y, y, 4);
		mpf_ui_sub(y, 1, y);
		mpf_sqrt(y, y);
		mpf_sqrt(y, y);
		mpf_add_ui(temp, y, 1);
		mpf_ui_sub(y, 1, y);
		mpf_div(yn, y, temp);
		
		// an = a*(1+yn)^4-2^(2*i+3)*yn*(1+yn+yn^2)		
		// Parte direita - 2^(2*i+3)*yn*(1+yn+yn^2)
		mpf_add_ui(temp, yn, 1);
		mpf_pow_ui(an, yn, 2);
		mpf_add(an, temp, an);
		mpf_mul(an, yn, an);
		mpf_mul_ui(an, an, twoexp);

		// Parte esquerda - a*(1+yn)^4
		mpf_pow_ui(temp, temp, 4);
		mpf_mul(temp, a, temp);
		
		// Parte esquerda - direita
		mpf_sub(an, temp, an);
		
		// Calcula pi anterior, pois o alg converge para 1/pi
		mpf_ui_div(temp, 1, a);
		// Atualiza anteriores
		mpf_set(a, an);
		mpf_set(y, yn);
		// Atualiza coef
		twoexp *= 4;
		
		// Calcula pi novo
		mpf_ui_div(an, 1, an);
		
		// Se nao mudou desde a ultima iteracao, escreve no arquivo e sai do loop
		if(mpf_cmp(temp, an) >= 0) {
			size_t writtenBytes = mpf_out_str(file, 10, 0, an);
			break;
		}		
	}
	// Limpeza
	close(file);
	mpf_clear(a);
	mpf_clear(y);
	mpf_clear(an);
	mpf_clear(yn);
	mpf_clear(temp);
	
	// Calculo do tempo decorrido
	gettimeofday(&final, NULL);
	time = final.tv_sec - inicio.tv_sec;
	printf("Tempo: %i segundos.\n", time);
	return 0;
}
